Out of the stable: Social disruption and concurrent shifts in the feral mare (Equus caballus) fecal microbiota

Abstract The disruption of animals' symbiotic bacterial communities (their microbiota) has been associated with myriad factors including changes to the diet, hormone levels, and various stressors. The maintenance of healthy bacterial communities may be especially challenging for social species as their microbiotas are also affected by group membership, social relationships, microbial transfer between individuals, and social stressors such as increased competition and rank maintenance. We investigated the effects of increased social instability, as determined by the number of group changes made by females, on the microbiota in free‐living, feral horses (Equus caballus) on Shackleford Banks, a barrier island off the North Carolina coast. Females leaving their groups to join new ones had fecal microbial communities that were similarly diverse but compositionally different than those of females that did not change groups. Changing groups was also associated with the increased abundance of a several bacterial genera and families. These changes may be significant as horses are heavily dependent upon their microbial communities for nutrient absorption. Though we cannot identify the particular mechanism(s) driving these changes, to the best of our knowledge, ours is the first study to demonstrate an association between acute social perturbations and the microbiota in a free‐ranging mammal.

diet, kinship, and shared environments), strongly implicating direct physical contact among social partners in the transmission of gut microbial species (Tung et al., 2015). Similarly, in semi-feral Welsh ponies (Equus ferus caballus), microbiota communities varied with spatial structuring and even with specific social interactions within groups (including those between males and females and mothers and offspring; Antwis et al., 2018).
In addition, various stressors including illness, environmental stressors, and physical trauma (experienced by virtually all species at one time or another) have been linked to microbial disruption in several taxa (Zaneveld et al., 2017). Factors as disparate as foraging in degraded habitats (howler monkeys, Alouatta pigra), pathogen infection (bullfrogs, Lithobates catesbeianus), and increases in temperature (corals, Anthozoa spp.) can all affect species' microbiota (Amato et al., 2013;Zaneveld et al., 2016). Social species experience similar stressors, in addition to those unique to group-living, including increased competition (for both food/water resources and mates), the maintenance of one's rank within and/or between groups, and group "politics" (Ward & Webster, 2016), all of which can potentially affect their microbiotas. Links between social stressors and the microbiota have been well-established in the laboratory. Mice (Mus musculus spp.) subject to prolonged aggression and infant rhesus monkeys (Macaca mulatta) subject to maternal separation show alterations to their microbiota (Bailey et al., 2010;Bailey & Coe, 1999). In the barn, domestic horses (Equus spp.) subject to transport or weaning and those with undifferentiated colitis (inflammation of the caecum and colon) all show changes to their microbiota (Costa et al., 2012;Mach et al., 2017;Schoster et al., 2016). Such alterations, often including decreased microbial diversity and stability, may reduce animals' ability to fight off infection (Bailey, 2012), leaving them more susceptible to a host of enteric pathogens (Bailey et al., 2010, Bailey & Coe, 1999. Similar links have been established in the wild. Both red squirrels (Tamiasciurus hudsonicus) experiencing trapping stress and yellow-legged gulls (Larus michahellis) with experimentally elevated corticosterone levels show decreased microbiota diversity.
Corals show myriad changes with simulated overfishing, nutrient pollution, and increased temperature (Noguera et al., 2018;Stothart et al., 2016;Zaneveld et al., 2016). Similarly, howler monkeys living in degraded habitats exhibit less diverse microbial communities; though these changes were attributed largely to dietary differences, the authors concede that increases in stress hormones might have played an important role (Amato et al., 2013).
Here, we investigate whether increased social instability is associated with changes in the microbiota in a population of free-living, feral horses (Equus caballus) living on Shackelford Banks, North Carolina, USA. Feral horses organize themselves into groups (bands) which contain the male(s)-or stallion(s), his female(s)-or mare(s), and their offspring (Feist & McCullough, 1976;Rubenstein, 1981Rubenstein, , 1986. These groups are normally long-lasting with few changes to group composition save for the dispersal of offspring (Klingel, 1975). Mares will leave their males/groups to join new ones primarily in response to a band stallion's attributes, including age, rank, length of tenure, and to a lesser extent, his propensity to engage in male-male aggression (Rubenstein, 1994;Rutberg, 1990), but at a fairly low rate (0%-17% of females; Berger, 1977, Feist & McCullough, 1976, Rubenstein, 1981. This stability of group membership is important to the overall well-being of all group members; decreases in stability have been associated with decreased body condition and reproductive success, and increased parasite load and offspring mortality (Kaseda et al., 1995;Linklater et al., 1999) all of which could affect animals' microbiotas (Comizzoli et al., 2021;Zaneveld et al., 2017). Moreover, the individual changing groups will experience a new home range and potentially changes in forage and water access, exposure to the microbial communities of their new bandmates, and a reassessment/reestablishment of their rank; again, all of which can affect their microbiota (Antwis et al., 2018;Zaneveld et al., 2017).
Decreases to the social stability typical of feral horses have been associated with immunocontraceptive treatment (Madosky et al., 2010;Nuñez et al., 2009Nuñez et al., , 2017Ransom et al., 2010). From 2000 to 2009, the feral horse population of Shackelford Banks, North Carolina was managed with the immunocontraceptive vaccine, porcine zona pellucida (PZP). Briefly, the vaccine induces an immune response which stimulates the production of antibodies that bind sperm receptors on the egg's surface, thereby blocking fertilization while still allowing ovulation and the associated cycles of sex hormones (Sacco, 1977).
Typically, anti-PZP titers only remain above control levels for approximately 10 months (Willis, 1994), though one PZP treatment can effect subfertility for up to three breeding seasons (Turner Jr. et al., 2007).
Mares treated with PZP can be 10 times as likely to switch bands, leaving their male to join another , thereby decreasing band stability and effecting some level of social disruption.
At the time of this study (6-7 years after the contraception program was largely suspended), Shackelford mares that previously received PZP continued to change groups more often compared to mares that had never been treated (Jones & Nuñez, 2019;Nuñez et al., 2017).
Evidence suggests that these prolonged effects are due to continued subfertility and not any lingering effect of the treatment itself (Nuñez et al., 2017). This increased variability in group changing behavior afforded us the unique opportunity to better discern the potential links between social disruption and microbiota in the wild.
In social species, significant changes in the behavior of individuals or subgroups of individuals often affects that of their close associates (Ward & Webster, 2016). Stallions respond to female group changing behavior aggressively with chases, and sometimes kicks and bites to the females as they attempt to bring them back to the group (Madosky, 2011).
Moreover, females engaging in this behavior often find themselves near highly escalated male-male conflicts (Jones & Nuñez, 2019). In addition, females that manage to change groups are often subject to aggression from resident females (C. M. V. Nuñez, unpublished data). Perhaps unsurprisingly, group changing behavior is associated with an endocrine stress response (Nuñez et al., 2014) which is reliably assessed via fecal cortisol metabolite (FCM) concentrations (Mostl & Palme, 2002;Nuñez et al., 2014;Wasser et al., 2000). Mares engaging in group switching exhibited higher FCM concentrations than those that did not. Moreover, mares engaging in the behavior more often exhibited higher FCM concentrations than those doing so less often. In some cases, FCM concentrations remained elevated for at least 2 weeks post-behavior (on average) indicating at least some lasting effect of the behavior (and/or its associated stressors) on mare stress physiology.
We sought to determine whether group switching behavior, that is, increased social instability, was correlated with changes in the mare's gut microbiota. Using fecal samples gathered in the field (2015 and 2016), we examined the microbiota communities of mares that did and did not change groups to determine links between a relevant social behavior (and potential stressor) and the diversity and composition of the microbiota. We predicted that mares changing groups would exhibit decreased diversity to their microbiota communities when compared to mares that did not change groups (Bailey et al., 2011). We also expected that mares changing groups would have compositionally different (shifted) fecal bacterial communities and would host bacterial communities that were more variable than those that did not (Zaneveld et al., 2017).

| Study area
We conducted this study on Shackleford Banks, a barrier island located approximately 3 km off the North Carolina coast, USA. The island is 15 km long and varies between 0.5 and 3 km in width.
There are approximately 14 dominant grass species that vary in distribution along the length of the island (Rubenstein, 1981)  The island comprises three distinct regions that also vary in habitat visibility (Jones & Nuñez, 2019;Rubenstein, 1981). The western portion of the island contains two primary water sources and is dominated by high dunes and dense brush that limit visibility; animals in this area share resources and exhibit overlapping home ranges.
The eastern portion is flat and open, with high visibility; historically, horses in this region were territorial (Rubenstein, 1981), though more recent decreases in available water have necessitated significant home range overlap (Jones & Nuñez, 2019;Marr, 1996). The central portion's terrain is similar to that of the east, facilitating high visibility, though water sources are distributed more evenly, allowing horses to maintain more bounded territories (Jones & Nuñez, 2019, Rubenstein, 1981. Approximately 29%, 32%, and 38% of the animals are located in the western, eastern, and central portions of the island, respectively. The horse population on Shackleford Banks has been comanaged by the National Park Service (NPS) and the Foundation for Shackleford Horses since 1996.

| Study subjects
Shackleford horses are typical of feral equids. They form coherent bands of one or sometimes two or three stallion(s) with one to several mare(s) and their offspring (Rubenstein, 1981). Although multimale bands can constitute a large portion of feral horse populations, (e.g., 33% in the Kaimanawa horses; , they occur less frequently on Shackleford Banks, accounting for only 10% and 5% of bands in the 2 years of this study (2015 and 2016, respectively). In addition, differences in the degree of stallion territoriality have been associated with the island's varying ecology (Rubenstein, 1981). During the time of this study, bands in the western portion of the island moved freely within overlapping home ranges, while bands in the central and eastern regions were more likely to defend territories (Jones & Nuñez, 2019). Regardless, bands living in all island regions remained spatially distinct from one another and band membership was easily determined (Feist & McCullough, 1976;Rubenstein, 1981Rubenstein, , 1986. Interactions between bands typically in- Historically, the bands on Shackleford Banks were stable with most changes to group membership involving the dispersal of immature individuals (Nuñez, 2000;Rubenstein, 1981). Females did leave their males/groups to join others, primarily in response to band stallion quality (Rubenstein, 1994;Rutberg, 1990), but at a low rate (10.8% of females; Rubenstein, 1981). Males sometimes fought to acquire mares from other groups, but stallions almost always retained their mares (Rubenstein, 1981) as is typical in other feral horse populations (Feist & McCullough, 1976). More recently, mares treated with PZP immunocontraception have changed groups more often, making approximately 10 times more group changes than untreated mares and joining (at least temporarily) five times as many groups (Madosky et al., 2010;Nuñez et al., 2009). At the time of this study, 6-7 years after suspension of contraception management, previously treated mares continued to change groups more often than mares that had never been treated (Jones & Nuñez, 2019;Nuñez et al., 2017), providing ample variation in this behavior for comparison. To control for any residual effects of PZP treatment to mare physiology and/or behavior (Nuñez et al., 2017), we considered only mares that had been contracepted at least once between 2000 and 2009 (S. Stuska, personal communication, May 2016); these mares received an average of 4.25 ± 0.28 total treatments (range = 1-7).

| PZP contraception
At the time of this study, PZP treatment had been suspended for approximately 98% of females: only two mares (one mare in 2015; that mare and one additional mare in 2016) actively received PZP treatment (see Table S1). These animals were contracepted primarily due to concerns that they would not survive pregnancy (S. Stuska, personal communication, May 2016).

| Animal welfare
All sampling was conducted in accordance with National Research  M. Fatkah and R. Schwartzbeck, 2016). All observers were trained by C.M.V. Nuñez. Comparable to previous studies (Nuñez et al., , 2017Ransom et al., 2010), data were collected across 8 weeks in 2015 (late June-mid-August) and 8.5 weeks in 2016 (mid-June-early August), totaling over 294 h of behavioral observation (192 h, 2015; 102.85 h, 2016), and averaging 3.39 ± 0.20 h per mare (range = 1.23-10.37 h). Horses were identified individually by color, sex, physical condition, and other distinguishing markings, including scars and freeze brands. Ages of the horses are known from long-term records (Nuñez, 2000;) and from NPS data (S. Stuska, unpublished data).

| Behavioral and demographic sampling
We located each study group once every 10 days and once every 3.5 days (on average) in 2015 and 2016, respectively. We recorded each group's GPS location and composition, noting the presence or absence of individual mares. We observed most of the reproductive mares (aged 4 years and older) present during the study, totaling 61 mares, representing 100% and 95% of the potential study population in 2015 and 2016, respectively. Fifty-four mares were observed in both seasons; the remaining mares were observed in either 2015 (n = 4) or 2016 (n = 3).
Only previously contracepted mares that produced fecal samples (n = 30) were considered for these analyses. Ten of the 30 mares that contributed samples changed groups during the study: two mares changed groups in both 2015 and 2016; three mares changed groups in 2015 only; five mares changed groups in 2016 only (see Table S2, Figure S1).
Mare transfer activity was rarely witnessed directly (2015, n = 1; 2016, n = 1); therefore, mare absence from a band was an important metric with which we measured the number of transfers between groups (Jones & Nuñez, 2019;Nuñez et al., 2009Nuñez et al., , 2017. We remained with each group for at least 30 min to ensure that individuals recorded as absent were not actually nearby, but out of our sight. Transfer behavior was confirmed by the mares' presence in new bands. In 2015, confirmation was made within 2.50 ± 0.35 days on average (range = 2-3); in 2016, confirmation was made within 13.58 ± 2.49 days on average (range = 1.5-33).

| Body condition
We assessed mare condition via rump scoring as physical condition has been associated with the microbial communities of horses (Garber et al., 2020). We examined the curvature of the line between the tailbone and the point of the hip to determine an individual's rump score. Scores were based on a scale from 1 to 5; a score of 1 being the poorest (Pollock, 1980). Females were scored an average of 1.33 ± 0.14 times (range = 1-2) and 2.08 ± 0.21 times (range = 1-4) in 2015 and 2016, respectively. For analysis, we scaled and centered scores at zero to aid interpretation of results (i.e., the resulting model intercepts can be interpreted as representing mares with an average rump score). For example, the parameter estimate for the intercept of a model of operational taxonomic unit (OTU) richness constructed (using unscaled average rump scores as a predictor variable) reflects the alpha diversity of a horse with a low rump score. Conversely, the intercept of a model constructed using the scaled rump score is representative of horses with an average rump score, a more biologically intuitive approach to model interpretation.

| Microbiota analysis
Fecal samples were sent to Argonne National Laboratories for DNA extraction, library preparation, and sequencing. DNA was extracted using the Qiagen DNeasy PowerSoil extraction kit (Qiagen Inc.). The V4 region of the 16S rRNA gene was amplified using the 515F and 806R primers (Caporaso et al., 2012), and barcoded using a primer set adapted for the Illumina MiSeq. Sequence data were generated using paired-end sequencing for 300 cycles using Illumina MiSeq v2 chemistry. We sequenced 2 × 150 bp amplicons on a MiSeq run, and samples from both years were sequenced in the same run.
Sequences were trimmed to 150 bp and assigned to suboperational taxonomic units with Deblur (Amir et al., 2017). We assigned taxonomy using the Silva 138 99% OTUs (515F/806R) classifier (Bokulich et al., 2018). A representative sequence for each OTU was aligned with mafft (Katoh et al., 2002) and used to construct a phylogeny with qiime phylogeny (Price et al., 2010) Note: No metric of alpha diversity (observed richness, Shannon-Weaver evenness, or Faith's phylogenetic diversity) differed between mares with respect to group change status; however, regional and temporal differences were detected.
(R Core Team, 2022). Prior to analysis of alpha diversity, we subsampled samples to 5706 reads (the number of reads present in the sample with the fewest reads, see Figure S2). Using the subsampled dataset, we calculated three alpha diversity metrics: richness (the number of observed OTUs in each sample), Shannon-Weaver evenness (a measure of the distribution of the relative abundance of each OTU in each sample), and Faith's phylogenetic diversity, which accounts for both the number of OTUs and their phylogenetic relatedness by summing the branch lengths among bacterial taxa represented in an individual's microbial community (Kembel et al., 2010). In our initial maximal models, we tested the effects of group changing behavior, study year, island region, female body condition (average rump score), and one interaction-the interaction between group changing behavior and body condition-on each alpha diversity metric. We tested for this interaction to address the possibility that changes in the microbiota associated with group changing behavior were influenced by mares' physical condition (Garber et al., 2020). All models included a random effect for horse ID to account for nonindependence introduced in situations where a single horse contributed multiple samples. Models were constructed to test all combinations of explanatory variables and were compared with Akaike's Information Criterion adjusted for small sample sizes (AIC c ; Burnham & Anderson, 2002). We report predictor variables retained in all models with ∆AIC c <2 (Table 1), and parameter estimates for the top models (∆AIC c = 0; Table 2).

Betadispersion
Betadispersion is a metric of variability in community composition, in this case, the dissimilarity or distance between an individual mare's microbial community composition and that of its group centroid. Betadispersion values can be used to assess the homogeneity of variance among sample groups (e.g., whether the compositional variance of the microbiota of mares that changed groups was equal to that of mares that did not change groups). We computed betadispersion for Bray-Curtis, UniFrac, and weighted UniFrac metrics using the "betadisper" function in the vegan package (Oksanen et al., 2018). Bray-Curtis dissimilarity reflects how similar (or dissimilar) samples or groups are based on the composition (identity and relative abundance, but not phylogenetic relatedness) of microbial taxa shared between samples. Higher scores in Bray-Curtis space indicate greater dissimilarity (two communities with a score closer to 1 share a smaller proportion of their species than two communities with a score closer to 0). Similar to Bray-Curtis dissimilarities, UniFrac distances account for presence and absence of taxa within communities accounting for phylogenetic relatedness among community members. If two communities differ in Bray-Curtis space but not UniFrac space, it may be because the taxa that differ between communities are closely related, having a high degree of shared branch lengths within the communities. Weighted UniFrac analysis provides a metric of communities' phylogenetic relatedness while accounting for the relative abundance of bacterial taxa, affording a comprehensive survey of the microbiota which is less sensitive to the presence of rare species than presence/absence-based metrics (Lozupone et al., 2007). Weighted UniFrac distances are quantitative representations of phylogenetic distances (branch lengths) and compositional differences (differences in relative abundances), between communities, with higher values indicating greater distances between communities (Lozupone et al., 2007).
We analyzed differences in betadispersion estimates using the "permutest" function in the vegan package and included a permutation structure to account for nonindependence of samples collected from the same mare. Betadispersion estimates were calculated with respect to median estimates for groups of mares that did or did not change groups.

| PERMANOVA
We compared microbial communities between mares that did and did not change groups (controlling for study year, island region, and female body condition) using Bray-Curtis dissimilarities, UniFrac distances, and weighted UniFrac distances calculated using phyloseq, and visualized our comparisons with nonmetric dimensional scaling (NMDS) ordinations in ggplot2 (Wickham, 2016). We conducted permutational analysis of variance (PERMANOVA) with the "adonis2" function (by = "margin") in the vegan package (Oksanen et al., 2018) to determine whether switching groups had any influence on mares' microbial communities. Our PERMANOVA also included variables to account for microbiota variation due to study year, island region, female body condition, and included a permutation design structure to account for nonindependence of samples collected from the same mare.

| Differentially abundant taxa
Finally, we analyzed whether any bacterial families or genera among the OTUs identified were differentially abundant in horse feces in association with several grouping variables including group changing behavior, study year, island region, and female body condition (when appropriate, see below). Estimating differential abundance among taxa in microbiota datasets is notoriously challenging because of the compositional nature of 16S rRNA data, and because differences in sampling fraction, the ratio of absolute abundance of a taxon in a single sample, versus its absolute abundance in the population of samples, can introduce bias and false positives (Lin & Peddada, 2020). We used two methods to assess differentially abundant taxa to accommodate the limitations inherent to each of the differential abundance estimation methods and our study's design. First, we used the more conservative ANCOM-BC method (Nearing et al., 2022) using the "ancombc" function from the R package ANCOMBC (Lin & Peddada, 2020) which accounts for both compositionality within samples, and differences in sampling fraction among samples. Using this method, we did not assess differentially abundant microbial taxa on the basis of horse condition because organizing this continuous metric into biologically relevant groups would have yielded unacceptably small sample sizes per group. Prior to testing for differential abundance with ANCOM-BC, we removed taxa observed in fewer than 10% of samples to avoid any undue influence of rare taxa in our analysis.
We then analyzed differential abundance on raw abundance data.
Next, we used the MaAsLin2 method from the Maaslin2 R package (Mallick et al., 2021) which other researchers have noted to be less conservative in comparison with ANCOM-BC, but highly consistent (Nearing et al., 2022), to allow for the inclusion of continuous predictor variables and random effects. Using this method, we were able to assess differentially abundant microbial taxa based on mare body condition. Additionally, we were able to include a random effect to account for instances in which multiple samples were collected from the same mares. We tested for differentially abundant taxa at the genus and family levels, again filtering out taxa that were present in fewer than 10% of samples, and whose relative abundances were <0.01%. Differential abundance analysis in MaAsLin2 was performed on raw (non-rarefied) count data, and data were transformed with the default log transformation.

| Alpha diversity
Group-changing behavior was not significantly associated with fecal microbiota richness, evenness, or phylogenetic diversity in mares (p obs = .333, p even = NA, p PPD = .257; Figure 1). Observed richness was best described by one model; Shannon-Weaver evenness was best described by an intercept-only model; and Faith's phylogenetic diversity was described by two models with ∆AIC c < 2. We report the parameter estimates, standard errors, and p-values of the top models (∆AIC c = 0) in Table 2.
Alpha diversity (calculated using richness and phylogenetic diversity metrics) was significantly lower in 2016 than in 2015 (p obs < .0001, p PD < .001). Additionally, samples collected from horses in the western region of the island tended to be less rich than samples from horses in the eastern and central regions of the island (p obs = .09). Phylogenetic diversity of western samples was also somewhat reduced in comparison to samples from other regions (p PD = .03, .11). No variation in community evenness was observed across years, regions, body condition, or group changing behavior status.

| Beta diversity
Betadispersion estimates for mares that changed groups were higher in both UniFrac and Bray-Curtis space (p = .002), but this pattern was not repeated with weighted UniFrac distances (Table 3, Figure 2).
In Bray-Curtis, UniFrac, and weighted UniFrac spaces, mares' fecal microbial communities were significantly different from each other based on group changing behavior even after controlling for the effects of study year, island region, female body condition, and repeated measurements of several individual mares (p all < .05; Table 4, Figure 3). Using all metrics of beta diversity, fecal bacterial communities also varied by year (p < .01; Table 4). Additionally, there was a significant association between female body condition and community composition in weighted UniFrac space (p = .008; Table 4).

| Differentially abundant taxa
Prior to running differential abundance analyses with ANCOM-BC

| Differentially abundant bacterial families
No bacterial families were shown to be differentially abundant in association with group changing behavior using ANCOM-BC, while F I G U R E 1 Group changing status of Shackleford Banks, NC mares was not associated with microbiota alpha diversity in terms of its richness (a), evenness (b), or phylogenetic diversity (c). Each point represents a single fecal sample; some mares produced multiple samples. Median alpha diversity for each metric is shown with a solid line, while upper and lower dashed lines represent the 75th and 25th quantiles, respectively.  (Table 5).

| Differentially abundant bacterial genera
Differential abundance analysis with ANCOM-BC identified four bacterial genera: an uncultured genus in the family Marinifilaceae, Caldicoprobacter (order Clostridiales, phylum Firmicutes), the NK4A136 group in the genus Lachnospiraceae, and an unnamed genus (p-1088-a5_gut_group) in family Pirellulaceae, that were significantly more abundant and one genus (Blautia) that was differentially less abundant in the feces of horses that changed groups (Table 5). Using MaAsLin2, we identified two genera, Eubacterium and an uncultured genus in family Erysipelotrichaceae that were significantly more abundant in feces from mares that changed groups, and one genus (Blautia) for which abundance was lower in mares that changed groups (Table 5). Only the bacterial genus Blautia was identified by both methods of differential abundance analysis as varying with group changing behavior (Figure 4).

| Associations with body condition, year, and region
Neither method of analysis identified differential abundances of taxa at the family or genus level with mare body condition. In contrast, both methods identified differential abundances of several taxa (at both the family and genus level) with study year and region: we detected nine differentially abundant families and seven differentially abundant genera associated with study year, and three differentially abundant families and two differentially abundant genera associated with study region (Tables S3 and S4).

| DISCUSS ION
Here, we show that fecal microbiota composition was associated with group changing behavior (a proxy for social instability) in feral mares.
Contrary to our predictions, no metric of alpha diversity, including richness, evenness, and phylogenetic diversity, differed for mares that changed groups compared to mares that did not. However, the composition of the gut microbiota of mares that changed groups was distinct from those of mares that did not in Bray-Curtis, UniFrac, and weighted UniFrac spaces. These results indicate that there is a small (R 2 = .025, .038) but significant (p < .05) difference between the composition of mares' fecal microbial communities that do and do not change groups.
We observed differences between communities using PERMANOVA, which is intended to detect significant shifts in group centroids, indicative of among-group variation. However, this type of analysis can be confounded by large differences in within-group variation. To address this possibility, we also compared betadispersion estimates among groups of mares and found differences in Bray-Curtis and UniFrac spaces associated with group changing behavior. Compositional differences in Bray-Curtis and UniFrac TA B L E 3 Betadispersion, or within-group variation in beta diversity, differed between mares that did and did not change groups when compared in Bray-Curtis and unweighted UniFrac space but did not differ between groups in weighted UniFrac space.

F I G U R E 2 Betadispersion estimates differed between
Shackleford Banks, NC mares that changed social groups and those that did not in Bray-Curtis space (a) and unweighted UniFrac space (b) but not in weighted UniFrac space (c). Each point represents a single fecal sample; some mares produced multiple samples. Median alpha diversity for each metric is shown with a solid line, while upper and lower dashed lines represent the 75th and 25th quantiles, respectively.

(a) (b) (c)
spaces could be driven by the heterogeneity of variances within communities, by differences in community membership, or by phylogenetic relatedness to varying levels. If greater within-group variation among mares that changed bands drives the communitylevel differences detected with PERMANOVA, it would suggest that social instability may precipitate microbiota community instability.
However, our results in weighted UniFrac space provide more conclusive evidence that bacterial communities of mares that did and did not change groups are compositionally different. Specifically, in weighted UniFrac space, the difference we observed betweengroup centroids was not accompanied by marked differences in within-group variation between mares that did and did not change social groups. Rather, our analyses in weighted UniFrac space provide evidence of a more directional shift in the composition of the fecal microbiota community. Taken together, our results suggest that social perturbations (like group changing behavior) may not precipitate wholesale shifts in bacterial diversity as measured by richness or relatedness (alpha diversity). Instead, social perturbations may initiate more subtle shifts in community structure. These changes are driven by differing levels of intragroup variation (depending on the metric used to compare communities), as well as combinations of changes in the relative abundances of bacterial taxa, the identities TA B L E 4 PERMANOVA analysis.

F I G U R E 3
Group changing behavior by Shackleford Banks, NC mares was associated with a shift in microbial community composition, as described by Bray-Curtis dissimilarities (a) and unweighted (b) and weighted UniFrac (c) distances. Each point represents a single fecal sample; some mares produced multiple samples. Stress values refer to the stress of the ordination solution and printed p-values are derived from PERMANOVA analysis describing the contribution of group-change status to microbiota variation (reported in Table 4).

(a) (b) (c)
of taxa within communities, and their relatedness within communities (beta diversity). Although our study does not address the consequences of such a shift, we argue that detecting these shifts and considering their possible mechanistic origins may shed light on the cascading effects of social perturbation in group-living animals.
There are myriad mechanisms that could be driving these changes. For example, shifts in microbial community composition are often linked to changes in diet (Amato et al., 2013;Ley et al., 2006Ley et al., , 2008Turnbaugh et al., 2009;Wu et al., 2011;Yildirim et al., 2010). In our study, dietary changes (including changes to both forage and water sources) and/or increased microbial transfer between newly encountered individuals (concomitant with changes in home range and/or region occurring after successful group switching behavior) likely played an important role in the changes to mares' fecal microbial communities (Stothart et al., 2021). Feral horse home ranges and group composition are typically stable (Berger, 1977;Feist & McCullough, 1976;Klingel, 1975;Rubenstein, 1981), resulting in more consistent diets and environmental microbiota, thereby driving microbiota similarity among group members (Antwis et al., 2018;Stothart et al., 2021). Increased dispersal between communities is thought to stabilize populations (Crowley, 1981); however, mares that change groups repeatedly may not remain in groups long enough to develop microbial communities that are representative of any one location (region or home range), instead experiencing increased exposure to bacterial communities from new environments and new group members. Work with semiferal Welsh ponies has demonstrated that microbiota communities can vary with spatial structuring and even with specific social interactions within bands (Antwis et al., 2018). Similarly, mares that change groups likely experience transitional fecal microbial communities from one band to the other, depending on how long they remain with each band; though again, these communities may not be representative of any one location. Such "transitional" microbiota communities have been demonstrated in dispersing male baboons: immigrant males that had lived in their current social group for a shorter period had microbiota that were less similar to other long-term group residents than did males with longer group residency times (Grieneisen et al., 2017). Such "transitional" microbial communities may leave mares that change groups (and particularly mares engaging in this behavior more often) more susceptible to environmental shifts; whether they constitute alterations to the surrounding environment (due to drought, hurricanes, increased/ decreased temperatures, etc.) or to the hosts' own physiology (due TA B L E 5 Differentially abundant bacterial families and genera detected in mare fecal samples from mares that changed groups versus those that did not, organized by method used for differential abundance analysis, and level of taxonomic identity (family or genus) and from most differentially absent (detected less than expected by chance) to differentially abundant (detected more than expected by chance).  (Table 6). Mean precipitation was slightly higher in 2015 (Table 6), though it is unlikely that this relatively small difference alone would have accounted for the stark differences in mare microbiota communities across years. While these conditions could indicate similar forage availability across years, shifts in diet cannot be ruled out. That said, even subtle differences in sociality have been shown to affect feral horse microbial community composition.
As previously described, studies with semi-feral Welsh ponies have demonstrated more similar fecal microbial communities among bands and even among more closely associated individuals within bands, indicating the importance of bacterial dispersal among close associates (Antwis et al., 2018). Work with the feral horses of Sable Island, Nova Scotia, demonstrated similar effects and revealed the importance of parental status (among mares) and ecological drift to feral horse microbiota communities (Stothart et al., 2021). It is possible that any one or several of these nonmutually exclusive effects could have influenced mares' microbiota communities differentially between years.
Still, even after accounting for year, island region, and mare body condition, the common act of changing groups was associated with changes to the microbial community, indicating that the social perturbation itself was associated with at least some of the changes to mares' microbiota communities. Group changing behavior can induce stress responses in Shackleford mares (Nuñez et al., 2014).
Mares changing groups are often chased, and sometimes kicked and/or bitten as resident stallions attempt to herd them back to the group (Madosky et al., 2010); they are subject to increased levels of reproductive behavior by both resident and new stallions ; and often find themselves near highly escalated malemale conflicts (Jones & Nuñez, 2019). Moreover, mares making successful changes to new bands are frequently subject to aggression from resident females (C. M. V. Nuñez, unpublished data). Though we cannot confirm that the differences exhibited here were due to mare stress physiology, chronic or prolonged stressors affect the microbial communities of diverse taxa in both the laboratory and in the wild.
Mice, rhesus monkeys, red squirrels, yellow-legged gulls, and corals all show changes to their microbiota with increased stress (Bailey et al., 2010(Bailey et al., , 2011Bailey & Coe, 1999;Noguera et al., 2018;Stothart et al., 2016;Zaneveld et al., 2016). Most relevant to this study, are the effects of transport, weaning, and colitis infection on the microbiota of domestic horses (Costa et al., 2012;Mach et al., 2017;Schoster et al., 2016). Myriad changes have been documented; however, decreases in the relative abundance of bacteria in order Clostridiales and increases in the abundance of bacteria in genus Bacteroides are notable (Costa et al., 2012, Mach et al., 2017, Schoster et al., 2016.
Interestingly, we found a pattern of consistent depletion of bacteria in the Blautia genus among mares that changed groups. The Blautia genus includes several bacteria that were previously classified in the genera Clostridium and Ruminococcus (Liu et al., 2021).
Bacteria in this genus are strictly anaerobic and some strains of this genus have been associated with probiotic activity and antiinflammatory activity (Liu et al., 2021). In captive horses, increased abundance of Blautia bacteria has been associated with an obesity phenotype (Biddle et al., 2018), while researchers studying human subjects in Japan noted a strong negative relationship between  Our previous work on Shackleford has shown that fewer contraceptive treatments, over longer periods of time can help to maintain more natural behavior and reproductive physiology (Nuñez et al., 2017) and reduce stress (Nuñez et al., 2014); given the present data, such management practice could also potentially maintain important symbiotic relationships in treated mares. Whether such relationships hold in other managed wildlife seems a worthwhile focus for future study. Conceptualization (lead); data curation (equal); formal analysis (supporting); funding acquisition (lead); investigation (supporting); methodology (equal); project administration (lead); resources (lead); supervision (lead); visualization (supporting); writing -original draft (equal); writing -review and editing (equal).

ACK N OWLED G M ENTS
We thank Dr. Sue Stuska with the National Park Service at the Cape Lookout National Seashore for her assistance in the field and for providing historical data on the horses. We also thank Maggie Kent, Kaylee Monroe, Micah Fatkah, and Rachel Schwartzbeck for their assistance in the field. Special thanks to Dr. James Adelman for his statistics advice and for his thoughtful comments on earlier versions of this article. We also thank the three anonymous reviewers who helped to improve this article. This work was supported by the

DATA AVA I L A B I L I T Y S TAT E M E N T
Genetic data: Sequence data are archived at the NCBI BioProject ID PRJNA962045 and can be accessed with the following link: https://www.ncbi.nlm.nih.gov/biopr oject/ 962045. Sample metadata: Metadata are available on the Open Science Framework data repository (Vaziri et al., 2022) at the following link: https://osf.io/rhz7a/.